'''
Created on May 21, 2011

@author: oabalbin
'''
import os
import subprocess
from optparse import OptionParser

def bed_tools_isec(bedtools_path, fileA, fileB, isecfile, local=False):
    '''
    Print intersection between file A and file B
    Print only unique entries infile A vs file B 
    '''
    
    bed_cmd = os.path.join(bedtools_path,"intersectBed")
    args = [bed_cmd, '-a', fileA, '-b', fileB, '-wa', '-u', '>', isecfile]    
    if local:
        outfile=open(isecfile,'w')
        args = [bed_cmd, '-a', fileA, '-b', fileB, '-wa', '-u']
        args = subprocess.Popen(args,stdout=outfile)
        outfile.close()
    
    return args

def get_header_lines(ifile,hederfile):
    '''
    '''
    output=open(hederfile,'w')
    p2 = subprocess.Popen(['grep','\#',ifile],stdout=output)
    output.close()
    return p2

def cat_files(afile, bfile, ofile):
    '''
    '''
    output=open(ofile,'w')
    p2 = subprocess.Popen(['cat',afile, bfile],stdout=output)
    output.close()
    return p2

    
def main_isec(bedtools_path, fileA, fileB, isecfile, 
              isecfile_header, hederfile,local):
    '''
    '''
    args = bed_tools_isec(bedtools_path, fileA, fileB, isecfile, local)
    p2 = get_header_lines(fileA,hederfile)
    p3 = cat_files(hederfile, isecfile, isecfile_header)
    
if __name__ == '__main__':
    optionparser = OptionParser("usage: %prog [options] ")
    # Specific options
    optionparser.add_option("-b", "--bed_tools_path", dest="bed_tools_path",
                            help="bed_tools_path")

    optionparser.add_option("-f", "--query_vcf_file", dest="query_vcf_file",
                            help="query_vcf_file")
    optionparser.add_option("-o", "--others_vcf_files",dest="others_vcf_files",
                            help="others_vcf_files. Pass vcf files without .gz extention even if the file .gz is there")
    
    (options, args) = optionparser.parse_args()
    print options
    
    local=True
    if not options.bed_tools_path:
        bedtools_path='/exds/users/oabalbin/sw/bedtools/current/bin/'
    else:
        bedtools_path=options.bed_tools_path
        
    path_to_exons_bed_file='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/refGene.B37.exons.pad50.bed'
    isecfile = options.query_vcf_file.replace('.vcf','_isec_SS50M.vcf')
    hederfile = options.query_vcf_file.replace('.vcf','_header.txt')
    isecfile_header=isecfile.replace('.vcf','h.vcf')
    
    if not options.others_vcf_files:
        fileB=path_to_exons_bed_file
    else:
        fileB=options.others_vcf_files
        
    main_isec(bedtools_path, options.query_vcf_file, fileB, isecfile, 
              isecfile_header, hederfile,local)
